Implementation of healthy eating policies across Europe

Introduction

Aim

In this paper, we use the NOURISHING framework and policy database to assess the implementation of healthy eating policies across Europe.

There is limited knowledge of implemented nutrition policies in Europe. The aim of this report is to showcase index construction by providing an overview of implemented policies in different European countries, and identify any similarities or clusters of policies between European countries. A simple regression displaying the index against obesity rates in European countries is included.

Limitations

• The WCRF data is not complete at this time.
• The current analysis does not assess quality or count of policies, only if a policy -AT ALL- is implemented in the given dimension at all. No judgement to the policy’s success is evaluated.
• Aggregation methods must be revisisted at a later time.

Findings

The findings indicate an overall low implementation of nutrition policies across European countries using the NOURISHING framework. A country may have multiple policies but an holistic approach to the policy landscape is missing.

For future index construction, the low policy coverage may pose a problem for an index using the benchmarks. This can affect aggregation, hence a scoreboard may be a better approach than an overall index. More investigation is necessary: correlation between aggregation levels, principal components analyses and sensitivity analyses.

Data

The data is manipulated in R prior to analyses. The following packages are utilised:

• Tidyverse is used to manipulate the data.
• COINr from the JRC is used to assemble and aggregate the data.
• Pheatmap is used for the cluster analyses.
• Reactable for making tables more user-friendly.
• RColorBrew for colourful tables.
• Countrycode for transforming country names into codes, and vice-versa.

WCRF policy database

Data is open source from http://policydatabase.wcrf.org.

policydb <- read_csv("C:\\Users\\JRMA\\OneDrive - Folkehelseinstituttet\\Dokumenter\\CO-CREATE\\policy-export10-Dec-2021.csv", show_col_types = FALSE)

colnames(policydb) <- c("DB_type", "PolicyArea", "SubPArea", "PolicyAction", "Country", "Topics", "BenchmarkID", "PolicyAreaID", "PolicyDim")

policydb <- policydb[-c(1,6)]

country_hierch2 <- policydb %>% count(Country, PolicyAreaID, BenchmarkID)

truef_spa <- country_hierch2 %>% mutate(
  PolicyActionImplemented = if_else(n > 0, TRUE, FALSE, missing=FALSE)
)

SPA_level_TF <- truef_spa %>%
  select("Country", "BenchmarkID","PolicyActionImplemented") %>% 
  select(c("Country", "BenchmarkID")) %>% distinct() %>%
  pivot_wider(names_from = "BenchmarkID", values_from = "BenchmarkID") %>%
  ungroup() %>%
  mutate(
    across(.cols = -Country, 
           ~if_else(!is.na(.), 1, 0)
    )
  )

AggHiearchy <- policydb %>% distinct(SubPArea, BenchmarkID, PolicyArea, PolicyAreaID)
PAreas <- policydb %>% distinct(PolicyArea, PolicyAreaID) #List of PAreas
SubPas <- policydb %>% distinct(SubPArea, BenchmarkID)

Descriptive statistics from the policy database

There are significant differences in number of policy actions collected per country:

In the following table it is evident that there are differences also by policy areas and sub-policy areas. The table is sort-able (e.g by PolicyArea).

So far, we know that there are differences in implemented of policy actions by both policy and sub-policy areas, as well as countries. This gives a clear indication that there will be differences in later stages of analyses.

Transform WCRF data to COINr package

The dataset is transformed to ensure consistency with the COINr framework. This is done to make a policy index where the policy dimensions are aggregated to a single number.

Indicator Data (IndData)

This an indicator of the indicators and units used the analyses. EEA countries are assigned to the EU, while EU countries not participating in the EU Comprehensive Scan outlined by WCRF is removed.

COINr_SPA <- SPA_level_TF

names(COINr_SPA)[names(COINr_SPA) == "Country"] <- "UnitName"


COINr_SPA <- COINr_SPA %>% mutate(UnitCode = countrycode(sourcevar = UnitName, origin="country.name", destination="iso3c"))

## Warning in countrycode_convert(sourcevar = sourcevar, origin = origin, destination = dest, : Some values were not matched unambiguously: Gulf Cooperation Council, Micronesia, Romania, Wallis and Futuna

## Warning in countrycode_convert(sourcevar = sourcevar, origin = origin, destination = dest, : Some strings were matched more than once, and therefore set to <NA> in the result: Romania, Wallis and Futuna,ROU,WLF

#Some values are missing: Gulf Cooperation Council, Micronesia, Romania, Wallis and Futuna

Missing <- COINr_SPA %>% subset(is.na(UnitCode)) %>% select(UnitName, UnitCode)
COINr_SPA <- na.omit(COINr_SPA)


isCoCREATE <- c("NOR", "NLD", "PRT", "POL", "GBR")

#Continue adding groups for countries 
COINr_SPA <- COINr_SPA %>% mutate(
  Group_Continent = countrycode(sourcevar = UnitName, origin="country.name", destination="continent"), 
  Group_EU = countrycode(sourcevar = UnitName, origin="country.name", destination="eu28"),
  Group_Region = countrycode(sourcevar = UnitName, origin="country.name", destination="un.regionsub.name"), 
  Group_COCREATE = case_when(UnitCode %in% unlist(isCoCREATE) ~ TRUE, TRUE ~ FALSE))

## Warning in countrycode_convert(sourcevar = sourcevar, origin = origin, destination = dest, : Some values were not matched unambiguously: Taiwan

COINr_SPA <- na.omit(COINr_SPA)

COINr_SPA <- COINr_SPA %>% mutate(Group_EU = case_when(
  UnitCode == "NOR" ~ "EU", 
  UnitCode == "ISL" ~ "EU", 
  UnitCode == "LIE" ~ "EU",
  UnitCode == "CHE" ~ "EU",
  UnitCode == "RUS" ~ "NA", #not part of scan this round
  TRUE ~ Group_EU)
)

COINr_SPA <- COINr_SPA %>% filter(Group_EU=="EU")


#only European countries for idnexation
COINr_SPA <- COINr_SPA %>% filter(Group_EU=="EU")

#Consider adding certain variables etc for further analysis (eurostat stuff)



IndData <- COINr_SPA

IndData <- IndData %>% mutate_if(is_double,as.integer)

IndData %>% reactable(compact = T)

Indicator metadata

The metadata for each indicators provides more details of single indicators and outlines the structure of the index.

IndMeta <- AggHiearchy

colnames(IndMeta) <- c("PolicyArea", "IndName", "IndCode", "Agg_PolicyDimension")
IndMeta$Direction <- 1 #All indicators got positive direction
IndMeta$IndWeight <- 1 #All indicators carry the same weight


#IndMeta$Agg1_Indicators <- IndMeta$IndCode ##Avventer denne. Blir dobbel uans? 


#IndMeta$Agg_PolicyDimension #Samles til en av de tre hovedkategoriene


#IndMeta$Agg_PolicyIndex #Alt samles

IndMeta$Agg_PolicyIndex <- "NOURISHING" 

IndMeta <- relocate(IndMeta, Agg_PolicyDimension, .before = Agg_PolicyIndex)
IndMeta %>% reactable(compact = T)

Aggregation metadata

This outlines how the index aggregates from indicator to full index (3 steps).

AggMeta <- PAreas

colnames(AggMeta) <- c("Name", "Code")

AggMeta$AgLevel <- 2

AggMeta$Weight <- 1 #Equal weight for all

AggMeta <- AggMeta %>% add_row(Name="Index", Code="NOURISHING", AgLevel=3, Weight=1)


AggMeta %>% reactable(compact = T)

Index construction

Assembling the data

Combine the IndData, IndMeta and AggMeta into one hierarchical list. This creates a “COIN” and will be the base we work from here on out.

In this dataset we have 67 indicators (benchmarks), 31 units (countries) and two aggregation levels (indicators aggregate to a policy area and then to the full index).

NPI <- assemble(IndData = IndData, IndMeta = IndMeta, AggMeta = AggMeta)



NPI

## --------------
## A COIN with...
## --------------
## Input:
##   Units: 32 (AUT, BEL, BGR, ...)
##   Indicators: 67 (G1, G2, G3, ...)
##   Denominators: 0 (none)
##   Groups: 4 (Group_Continent, Group_EU, Group_Region, ...)
## 
## Structure:
##   Level 1: 67 indicators (G1, G2, G3, ...) 
##   Level 2: 10 groups (G, H, I, ...) 
##   Level 3: 1 groups (NOURISHING) 
## 
## Data sets:
##   Raw (32 units)

Structure

The sunburst plot shows the hierarchical structure and effective weights for the simplified NOURISHING index.

Effective weights imply their relative contribution to the aggregated level, even with equal weighting applied. This is most apparent at indicator level: the N (level 2) has eight indicators contributing 0.0125 each, while in N(2), one indicator contributes 0.033. At sub-index level (N, O, U, R, I, S , H, I, N , G) they are contributing equally at 0.1 each.

Example of map

An example of countries implementing R1:

Descriptive statistics of raw data

The data indicate that some sub-policy areas are more covered than others (sort by “mean”).

The table indicate which areas have the ‘best’ coverage. This is I(2)1, N1 and N7 respectively. Multiple variables have low coverage (sort by MEAN), e. g. G5, H4 and R8.

Re-scaling data

Normalisation is not necessary if it is all on same scale, but moving from a smaller to larger scale will simplify presentatation of the results.

There are multiple ways to rescale data, but due to the binary nature of all indicators, a multiplication by 100 is sufficient. After this normalisation, the scores will either be 0 or 100 instead of 0-1.

NPI <- normalise(NPI, dset="Raw", ntype="custom", npara = list(custom = function(x) {x*100}))

Aggregation

A central step to create an index is the aggregation method. We do not want total compensability because being good in one area should not compensate for implementing less in another area.

In the first aggregation step, from indicator to policy area (e.g: from R1 to R) we use the arithmetic mean (simple averages).

The aggregation can be read left to right. The last column, NOURISHING, displays the total score for the index. Results are presented on the next page.

NPI <- aggregate(NPI, dset = "Normalised")

#A snapshot of the aggregation sequence
NPI$Data$Aggregated[(ncol(NPI$Data$Aggregated)-10): ncol(NPI$Data$Aggregated)] %>% roundDF(1) %>% reactable(compact = T)

Index results

The ranking of each country is in the plot and table below. Highlighted countries are members of the CO-CREATE project.

Barplot of scores

Ranking table

World map

Dimension scores by colour

The countries perform differently according to policy areas as indicated by a deeper, green colour:

Spidercharts

Examples using CO-CREATE countries. The overall mean for the index is (17.6475), but the implementation varies between the policy areas.

The spider charts are clickable. The first shows the overall mean, the second the mean of CO-CREATE countries.

Overall mean

iplotRadar(NPI, dset ="Aggregated", usel = isCoCREATE, aglev=2, addstat="mean") 

CO-CREATE mean

iplotRadar(NPI, dset ="Aggregated", usel = isCoCREATE, aglev=2, statgroup = "Group_COCREATE", addstat = "groupmean")

Cluster analyses

The vast number of indicators and dimensions may make it hard to see patterns. The following analyses, we try to make two clusters: one for benchmarks and one for policy areas (as columns), with countries as the row.

Policy area analysis

#Prepare data
Aggregates <- getResults(NPI, tab_type = "Aggregates") 

Aggregates <- Aggregates %>% mutate(
  Group_Continent = countrycode(sourcevar = UnitName, origin="country.name", destination="continent"), 
  Group_EU = countrycode(sourcevar = UnitName, origin="country.name", destination="eu28"),
  Group_Region = countrycode(sourcevar = UnitName, origin="country.name", destination="un.regionsub.name"), 
  Group_COCREATE = case_when(UnitCode %in% unlist(isCoCREATE) ~ TRUE, TRUE ~ FALSE))

Aggregates <- Aggregates %>% mutate(Group_EU = case_when(
  UnitName == "NOR" ~ "EU", 
  UnitName == "ISL" ~ "EU", 
  UnitName == "LIE" ~ "EU",
  UnitName == "CHE" ~ "EU",
  TRUE ~ Group_EU)
)

#Data without text
NumericalsOnly <- Aggregates

rownames(NumericalsOnly) <- sapply(Aggregates$UnitName,function(x) strsplit(as.character(x),split = "\\\\")[[1]][1])

NumericalsOnly <- NumericalsOnly %>% select(-Rank, -NOURISHING, -UnitName, -UnitCode, -Group_EU, -Group_COCREATE, -Group_Region, -Group_Continent)

#Make unique group for visualisation
Regions <- Aggregates %>% select(Group_Region)  

rownames(Regions) = rownames(NumericalsOnly) 

The results for each policy area is set to divide into five categories of countries and three categories of policy areas.

• UK is the clear ‘winner’ in all policy areas, but are weakest in the policy areas U and R. They are strongest in N.2, followed-by = and I.2.
• The country with the weakest implementation seems to be Bulgaria.
• CO-CREATE countries perform strongly with the exception of Poland and the Netherlands.
• Country region does not seem to influence implementation at first glance.
• It is evident that letter N is the most implemented of all policy areas (1st column)
• The policy areas G, S, H, U and R are the least implemented areas overall for the European countries.
• N.2. is not covere at all by over half the countries included in the analyses.

NumericalsOnly2 <- NumericalsOnly %>% roundDF(decimals=1)

NumericalsOnly2[NumericalsOnly2 < 0.01] <- NA
NumericalsOnly2[is.na(NumericalsOnly2)] <- ""

library("RColorBrewer")

## Warning: package 'RColorBrewer' was built under R version 4.1.1

pheatmap::pheatmap(NumericalsOnly, cluster_cols = T, main="Clustering of countries by policy area", cutree_rows=3, cutree_cols=3, annotation_row = Regions, display_numbers = NumericalsOnly2, color=colorRampPalette(c("snow3", "yellowgreen", "forestgreen"))(50))

Sub-policy area analysis

#Prepare data
AggregatesSubP <- getResults(NPI, tab_type = "Full")

#Data without text
NumericalsSuBPOnly <- AggregatesSubP

rownames(NumericalsSuBPOnly) <- sapply(NumericalsSuBPOnly$UnitName,function(x) strsplit(as.character(x),split = "\\\\")[[1]][1])

NumericalsSuBPOnly <- NumericalsSuBPOnly %>% select(-Rank, -NOURISHING, -UnitName, -UnitCode, -Group_EU, -Group_COCREATE, -Group_Region, -Group_Continent, -N, -O, -U, -R, -I, -S, -H, -"I.2.", -"N.2.", -G)

##ALTERNATIV
#Due to Pheatmap categories, I need to transpose the table to get the column annotation
trans_NumericalsSuBPOnly <- NumericalsSuBPOnly %>% t() %>% as.data.frame()

Dimensions <- AggHiearchy %>% select(BenchmarkID, PolicyAreaID) %>% as.data.frame()

Dimensions <- Dimensions %>% mutate(across("BenchmarkID", str_replace, "\\(", "."))

Dimensions <- Dimensions %>% mutate(across("BenchmarkID", str_replace, "\\)", "."))


#Add PolicyID and BenchmarkID
trans_NumericalsSuBPOnly <- left_join(trans_NumericalsSuBPOnly %>% mutate(BenchmarkID = rownames(trans_NumericalsSuBPOnly)), Dimensions, by = "BenchmarkID") 

trans_NumericalsSuBPOnly <- trans_NumericalsSuBPOnly %>% column_to_rownames("BenchmarkID")

SubPAz <- trans_NumericalsSuBPOnly %>% select(PolicyAreaID)

NumericalsSuBPOnly <- trans_NumericalsSuBPOnly %>% select(-PolicyAreaID) %>% t() %>% as.data.frame()

When using the sub-policy areas to explore differences, the country ranking changes.

Important! These variables are “Yes” or “No”: hence the colour is either red for Implemented, or blue for Not implemented.

pheatmap::pheatmap(NumericalsSuBPOnly, cluster_cols = T, cluster_rows = T, main="Clustering of countries by benchmarks", cutree_rows=4, cutree_cols=3, annotation_col = SubPAz, legend = F, color = colorRampPalette(c("snow3", "forestgreen"))(50))

Regression

There are multiple approaches to using the final index score. An important prerequisite is a sound construction of the index. This will take some time to ensure and may end up not being feasible.

The regression can be used for a specific policy area (e. g., Nutrition labelling or Restricting marketing), or the full index. Both approaches should be explored. Note, there is less variance in some areas than others that we need to investigate closer.

Get data on obesity

We can retrieve any data from Eurostat or other sources. Here, I use Obesity statistics for the full population in European countries.

library(stringr)

hlth_ehis_bm1e <- eurostat::get_eurostat("hlth_ehis_bm1e")

## Table hlth_ehis_bm1e cached at C:\Users\JRMA\AppData\Local\Temp\RtmpMdZAhN/eurostat/hlth_ehis_bm1e_date_code_FF.rds

BMI_cat <- hlth_ehis_bm1e %>% select(bmi) %>% distinct()

hlth_ehis_bm1e$year <- str_sub(hlth_ehis_bm1e$time,1,4)

BMI_euro <- hlth_ehis_bm1e %>% filter(sex=="T" & age=="TOTAL" & isced11=="TOTAL" & bmi=="BMI_GE30" & year=="2019") %>% select(geo, values)

BMI_euro <- BMI_euro %>% mutate(geo2= countrycode(sourcevar = geo, origin="eurostat", destination="iso3c")) %>% select(geo2, values)

## Warning in countrycode_convert(sourcevar = sourcevar, origin = origin, destination = dest, : Some values were not matched unambiguously: EU27_2020

colnames(BMI_euro) <- c("UnitCode", "Obesity rate")

BMI_euro %>% reactable::reactable()

Link index with data

After retriving the data from Eurostat, we combine the data from NOURISHING and the selected data source. Note, LIE has missing data. Data from different sources imported for GBR and CHE. For a final analysis, we should use HSBC or something else.

IndexScore_geo <- getResults(NPI, tab_type = "Summ") %>% select(UnitCode, NOURISHING)

IndexScore_geo <- left_join(IndexScore_geo, BMI_euro, by = c("UnitCode" ="UnitCode"))

#Add some missing data for this exercise. This data is not cross-checked!
IndexScore_geo <- IndexScore_geo %>% mutate("Obesity rate" = replace(`Obesity rate`, UnitCode=="GBR", 28))
IndexScore_geo <- IndexScore_geo %>% mutate("Obesity rate" = replace(`Obesity rate`, UnitCode=="CHE", 11.3))

IndexScore_geo %>% reactable::reactable(compact=T)

Regression line

Creating a fitted line we find that there is a positive relationship between higher implementation of nutrition policies and obesity rates.

No significant correlation between index and obesity rates identified (p=0.35).

regline <- IndexScore_geo %>% filter(!is.na(`Obesity rate`)) %>% lm(NOURISHING ~ `Obesity rate`,.) %>% fitted.values()

corell <- cor.test(IndexScore_geo$NOURISHING, IndexScore_geo$`Obesity rate`, method = "pearson")

corell

## 
##  Pearson's product-moment correlation
## 
## data:  IndexScore_geo$NOURISHING and IndexScore_geo$`Obesity rate`
## t = 0.94849, df = 29, p-value = 0.3507
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1927261  0.4972377
## sample estimates:
##       cor 
## 0.1734597

x <- IndexScore_geo$NOURISHING
y <- IndexScore_geo$`Obesity rate`

Regression plot

Create a scatterplot with NOURISHING score on the vertical axis, and obesity rate on the horizontal axis:

meanNourish <- mean(IndexScore_geo$NOURISHING)

plot <- IndexScore_geo %>% filter(!is.na(`Obesity rate`)) %>%
  plotly::plot_ly(x=~`Obesity rate`, y= ~NOURISHING, mode="markers", text =~UnitCode) %>%
  plotly::add_markers(y= ~NOURISHING) %>%
  plotly::add_trace(x=~`Obesity rate`, y=regline, mode="lines") %>%
  plotly::add_lines(y=meanNourish) %>%
  plotly::layout(showlegend=F) %>%
  plotly::add_text(textposition="top right")


plot

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter

Annex 1: Codes